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Abstract 

The dynamics of loop formation by linear polymer chains has been a topic of several theo- 
retical/experimental studies. Formation of loops and their opening are key processes in many 
important biological processes. Loop formation in flexible chains has been extensively studied by 
many groups. However, in the more realistic case of semiflexible polymers, not much results are 
available. In a recent study (K. P. Santo and K. L. Sebastian, Phys. Rev. E, 73, 031293 (2006)), we 
investigated opening dynamics of semiflexible loops in the short chain limit and presented results 
for opening rates as a function of the length of the chain. We presented an approximate model for 
a semiflexible polymer in the rod limit, based on a semiclassical expansion of the bending energy 
of the chain. The model provided an easy way to describe the dynamics. In this paper, using 
this model, we investigate the reverse process, i.e., the loop formation dynamics of a semiflexi- 
ble polymer chain by describing the process as a diffusion-controlled reaction. We make use of 
the "closure approximation" of Wilemski and Fixmann (G. Wilemski and M. Fixmann, J. Chem. 
Phys., 60, 878 (1974)), in which a sink function is used to represent the reaction. We perform 
a detailed multidimensional analysis of the problem and calculate closing times for a semiflexible 
chain. We show that for short chains, the loop formation time r decreases with the contour length 
of the polymer. But for longer chains, it increases with length obeying a power law and so it has a 
minimum at an intermediate length. In terms of dimensionless variables, the closing time is found 
to be given by r ~ L n exp (const. /L), where n = 4.5 — 6. The minimum loop formation time occurs 
at a length L m of about 2.2 — 2.4. These are, indeed, the results that are physically expected, but 
a multidimensional analysis leading to these results does not seem to exist in the literature so far. 
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I. INTRODUCTION 



In a previous paper [l], the dynamics of opening of a weak bond between the two ends 
of a semiflexible polymer chain was considered in detail. An approximate model for a stiff 
polymer ring in the rod limit, based on a "semiclassical" method was developed. This model, 
though approximate, was found to provide an easy approach to describe the dynamics of a 
worm-like polymer chain in the rod limit. In Ref. we used the model to analyze the 
dynamics of opening and to calculate the rates of opening as a function of length in the 
short chain limit. Here in this paper, we analyze the dynamics of loop formation. 

The closing dynamics of polymer chains has been studied extensively, being the key 
process in important biological functions, such as control of gene expression (2), [3], DNA 
replication 4| and protein folding. Experimental studies on loop formation involve monitor- 



ing the dynamics of DNA hairpins 



ling, lijxperimen 
irpins 



[9] and small peptides 



0, Q Q, 



Q, using 



fluorescence spectroscopic techniques. Several theoretical approaches are available for ana- 
yzing loop formation of a flexible chain. Using the formalism of Wilemski and Fixman (WL) 
14| for diffusion controlled reactions, the closing time r for a flexible chain was calculated by 
Doi and was found to vary as, r ~ L 2 . In another important approach, Szabo, Schulten 
and Schulten (SSS) [16| calculated the mean first passage time for closing for a gaussian 
chain and found r ~ L 3 / 2 . The two approaches have been analyzed by recent simulations 
[17I . \\\ . But real polymers such as DNA, RNA and proteins are not flexible and hence, 
it is more important to understand the closing dynamics of stiff chains. Unfortunately, in 
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24j are available in the 



this case, only simple, approximate approaches 
literature so far. Since worm-like chains are represented by differentiable curves, one has 
to incorporate the constraint |u(s)| = 1 and this has been a problem in dealing with semi- 
flexible polymers. Yamakawa and Stockmayer 25| and Shimada and Yamakawa 26j have 
calculated the static ring closure probabilities for worm-like chains and helical worm-like 
chains. According to their analysis, the ring closure probability for a worm-like chain has 
the form, G(0;L) = 896.32(/ p /L) 5 exp(-14.054/ p /L + 0.246L/Z P ), where l p is the persistence 
length of the chain. An approximate treatment that leads to the end-to-end probability 
distribution for semiflexible polymers has been given by Winkler et. al [2jJ and using their 



approach, the closing dynamics has been analyzed recently by Cherayil and Dua [27[ . They 



find that the closing time r ~ L u , where v is in the range 2.2 to 2.4. In an interesting 
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paper, Jun. et. al [28| showed that the closing time should decrease with length in the short 
chain limit and then increase with length for longer chains. Hence, the closing time has a 
minimum at an intermediate length. The reason for this behavior is that, for short chains, 
the bending energy contributes significantly to the activation energy for the process. Thus 
the activation energy ~ const. / L and therefore the closing time r ~ exp (const. /L). For 
longer chains, the free energy barrier for closing is due to the configurational entropy and 
hence, r obeys a power law. Jun et. al have followed an approximate one dimensional 
Kramers approach to reproduce this behavior and obtain the minimum of closing time at a 
length L m = 3.4/„, where l p is the persistence length of the chain. Monte Carlo simulations 
by Chen et.al lead to L m = 2.85/ p . See also the paper by Ranjith et al 3p| . 

In Ref. jl|, we analyzed the opening dynamics of a semiflexible polymer ring formed 
by a weak bond between the ends. We developed a model that describe the polymer near 
the ring configuration, using a semiclassical expansion of the bending energy of the chain. 
The model, though approximate, provided an easy way to analyze the dynamics. Using this 
model, we calculated the opening rates as function of the contour length of the chain. The 
formalism presented in Ref. 1|, took into account of the inextensibility constraint, |u(s)| = 1 
for semiflexible chains rigourously. The conformations of the chain can be mapped onto 
the paths of a Brownian particle on a unit sphere. We performed a semiclassical expansion 
about the most probable path assuming that the fluctuations about the most probable path 
are small. For the ring, we took most probable path to be the great circle on the sphere. 
This is again an approximation, as the minimum energy configuration for a semiflexible 
polymer loop does not correspond to the great circle. However, as described in Ref. [jl, it 
led to minimum energy values very close to exact results by Yamakawa and Stockmayer [25( 



and the approximation scheme by Kulic and Schiessel [311 ] . Once the ends of a semiflexible 
polymer are brought together, they can separate in any of the three directions in space. Our 
analysis showed that two of the three directions in space are unstable, while one direction is 
stable. If one considers the ring to be in the XY-plane, with its ends meeting on the Y axis, 
then the motion that leads to separation along the Y direction is stable, while the motions 
lead to separation along X or Z direction are unstable. The nature of instabilities along 
the X and Z directions are different. Hence, near the ring, the three directions in space are 
non-equivalent for a semiflexible polymer and are governed by different energetics (see Sec. 
IIII Bl) . One may also perform the expansion near the rod configuration by expanding about 



4 



the straight rod. On the unit sphere the straight rod corresponds to a point and unlike the 
great circle this is an exact minimum energy configuration (see Sec. IHIj) . 

In this paper, we present a detailed multidimensional analysis of the dynamics of loop 
formation in semiflexible chains. We make use of the approximation scheme developed in 
Ref. jj. Following Wilemski and Fixman [14j, the looping is described as a diffusion- 
controlled reaction. In the WF theory, the effect of the reaction is incorporated into the 
model using a sink function. In special cases, exact analytical results are possible for a 
delta function sink 
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34|. But for an arbitrary sink, and multidimensional dynamics, 



this is not possible. For such cases, WF suggested an approximation known as the "closure" 
approximation. In this, the diffusion limited life-time of the process is expressed in terms of a 
sink-sink correlation function and the essential step for finding the loop formation time is to 
calculate this sink-sink correlation function. For this, we need to know the time-dependent 
Green's function of the chain and the equilibrium probability distribution. We therefore 
derive the time-dependent multidimensional Green's function of the semiflexible polymer 
near the loop configuration by performing a normal mode analysis. This Green's function 
is then used to find the sink-sink correlation function for a Gaussian sink and the closing 
time. We find that the closing time r ~ (L/l p ) n exp(Al p /L), The exponent n = 4.5 - 6. r 
is found to be a minimum at a length L m ~ 2.2 — 2Al p which has to be compared with the 
value 3.4/p obtained in Ref. 2j]] and 2.85/ p of Ref. j^. We find L m to be weakly dependent 
on the range of the interaction between the ends. Thus, our analysis leads to results that 
are physically expected. It is worth mentioning that a multidimensional analysis leading 
to these results does not seem to exist in the literature so far. We also calculate the loop 
formation probability G(0; L) and find that our method leads to the correct behavior, i.e., 
G(0; L) ~ L~ 5 exp(— const. /L), thus showing that the procedure reproduces the previous 
results for this quantity [261 ]. 

The paper is organized as follows: In SecHT], we give a summary of the WF theory for 
diffusion controlled reactions and the "closure" approximation. In Sec. IHIl the semiclassical 
approximation scheme for bending energy of a semiflexible polymer is briefly outlined. The 
time dependent Green's function of the polymer is derived though a normal mode analysis 
near the loop configuration in Sec. IIVI The approximate probability distribution of the 
chain is given in Sec. IIVI In Sec. IVI A\ we calculate the sink-sink correlation function for a 
Gaussian sink and the closing time. In Sec. IVI Bl we give numerical results. Summary and 
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conclusions are given in Sec. IVHI 



II. THE "CLOSURE" APPROXIMATION 

In this section, we summarize the theory of diffusion-controlled intra-chain reactions of 
polymers developed by Wilemski and Fixman [l4| and their " closure" approximation for an 
arbitrary sink function. The dynamics of a single polymer chain in a viscous environment 
is governed by the diffusion equation, 

OP 

~St + Dp = w 

where D is the diffusion operator for the chain. If the chain is represented by N + 1 beads 
with position vectors represented by r = (ri, r 2 , ....rjy+i), then the general form of the 
diffusion operator is given by 

D = D J2 Vi.(V< + feT)- 1 ^). (2) 

i=i 

D = fcgT/£, is the diffusion coefficient of the segments and £ is the friction coefficient of the 
segments. Fj = VjZ7, where U is the potential energy of the chain. Eq. ([1]) may be solved 
to obtain the equilibrium distribution P eq of the chain, which is time independent. But if 
the chain has reactive ends, they can react and form a loop when they come sufficiently 
close and hence the probability distribution of an open chain will decay in time. In such a 
case, one may solve Eq. (TjQ) with appropriate boundary conditions. An alternate approach 
to the same problem is to introduce a sink function into the equation for P(r, t) as done by 
Wilemski and Fixman [14] . Then the reaction-diffusion equation that governs the dynamics 
of a polymer chain with reactive ends is 

^ + DP = -k r S(v)P, (3) 

where P is the distribution function of the open polymer chain and k r is the strength of 
the sink function. k r determines the rate at which the reaction occurs when the ends are 
sufficiently close. S is the sink function and is a function of rx,r 2 , ...rjv+i- Integrating Eq. 
(j3J) over all the coordinates r we get 

= -k r v(t), (4) 



where 



and 



v(t) 



J drS( 



r)P(r,t) 



(5) 



P s (t) = J drP{r,t) (6) 
is the survival probability. The function S can be any suitable function, but is usually taken 



to be a delta function or a gaussian. Eq. (j3J) can be solved exact 
for a delta function sink or a quadratic sink (see references 
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y only in one dimension 



341 ] , and the references 



therein). Therefore, WF introduced the assumption that P(r,t) may be approximated as 



P(r,t) = P eq (r)u{t), 



where 



v(t) 



(7) 



(8) 



■eq 



with 

v eq = j dvS{v)P eq {v). (9) 

This is referred to as the "closure" approximation. The average time of closing is the integral 
of the survival probability and is given by 



P s (t)dt = P s (0), 



(10) 



where P s {s) is the Laplace transform of P s {t). r is also expressed in terms of a sink-sink 
correlation function and in the diffusion-limited (k r — > oo) limit, it is given by 14j 



T 



V{t) 



eq 



1 dt. 



(11) 



T>(t) is the sink-sink correlation function defined by 

V(t) = J dr J dr>S(r)G (r,r';t)S(r')P eq (r>), (12) 

where Gq(t, r'; t) is the Green's function for the diffusive motion of the chain in the absence 
of the sink. Eq. (ITT]) was obtained by WF [14| . Note that V(oo) = v \ . To calculate 
T>(t), one needs to know Go(r,r';t) and -P e g( r ) an d these will be calculated in the following 
sections. We shall take S(r) to be a Gaussian, given by 



S(r) = S(R) = S x {R x )S y {R y )S z {R z ), 



(13) 



where R is the end to end vector for the chain and 



Si(Ri) = e-Wfrft/lVto rj), i = x,y or z. (14) 
T) is the width of the Gaussian sink. 

III. THE SEMI-CLASSICAL APPROXIMATION SCHEME FOR THE BENDING 
ENERGY 

In Ref.jl], we introduced an approximation scheme for the bending energy of a semiflexible 

polymer ring, which is based on a "semiclassical" expansion. In this section, we give a brief 

account of the approach. A semiflexible polymer is usually considered as a continuous, 

inextensible space curve represented by the position vector r(s), where s is the arc- length 

parameter. The bending energy of the chain is given by 

L 2 
k f ( d r(s^ ~ 



Em = 2 J W ) " S - (15) 



k is the bending rigidity. Since the curve is differentiable, one has the constraint, 

|u(s)| = l, (16) 

where u(s) = dr(s)/ds, the tangent vector at the point s. The partition function of the 
semiflexible polymer is the functional integral over the conformations represented by r(s), 

n f \ f -E bend [r(s)}\ 

Dr(s)ex H k B r )■ (17) 

This functional integral has to be performed with the constraint of Eq. f|T6l) . However, 
incorporating this constraint has been a problem in dealing with semiflexible polymers. In 
Ref. [ij, we wrote the partition function as an integral over u(s), 

Z = y D U ( g )exp(, k B T J' (18) 
and represented u(s) in angle coordinates 

u(s) = isin#(s) cos0(s) + j sin6 l (s) sin0(s) + kcos6 l (s). (19) 



S 



FIG. 1: The conformations of a semiflexible polymer may be mapped onto the paths of a Brownian 
particle on the surface of a unit sphere. The circular ring polymer with the tangent vectors at the 
ends joining smoothly corresponds to the great circle on the unit sphere. 



Since the magnitude of the tangent vector is one, the conformations of the semiflexible 
polymer can be mapped onto the trajectories of a Brownian particle over a unit sphere 
(FigUJ). The bending energy of the chain is then written in terms of the angles 9 and 4> as 

L 

E bend ='- fds\( + sin 2 0( S ) ( =rp. ) \ (20) 





and the partition function is written as a path integral in spherical polar coordinates 

Z*. = / DOW J i*W exp (- ^y^ ) • (21) 

This path integral has not been evaluated in a closed form. In Ref. we have used 
a semiclassical expansion of the bending energy, to evaluate the above partition function 
approximately. 



A. Bending energy of the loop: expansion about the great circle 

To perform a semiclassical expansion of the bending energy of the semiflexible polymer 
near the ring configuration, we take the most important path to be the great circle on the 

9 



unit sphere (Fig. [T]). The great circle corresponds to a ring with the tangents smoothly 
joined. However, the minimum energy configuration of a rod-like polymer whose ends are 
brought together to form a loop would not have its tangents joining smoothly and therefore 
does not correspond to a great circle 2^|. Hence, our approach is approximate but has the 
advantage that it provides an easy way to study the dynamics. On the other hand, if one 
is interested in covalent bond formation, in which directionality of the bond is important, 
then the great circle is the appropriate starting point. 

The position vector of the polymer may be found by inverting the definition of the tangent 
vector u(s) = dr(s)/ds, 



s 



r(s) = r cm - — / ds dsiu(s x ) + / dsiu(s x ), (22) 
L> Jo Jo Jo 

where r cm denotes the position vector of the center of mass of the ring polymer. The great 
circle is chosen to lie in the XY plane of a Cartesian coordinate system, with any point on 
it represented by the coordinates, [9(s),<p(s)] = [7r/2,2irs/L]. The position vector of the 
circular ring polymer that corresponds to the great circle may be found using Eq. (1221) and 
is given by 

T, I 9'7re\ / 9tt e \ 

(23) 



27TS\ . ( 2lTS 

i sin l — — — j cos 



L J \ L 

This curve represents one end of the polymer lying in the XY-plane starting at ^ on the 



2tt 

negative Y axis, going around the Z axis along a circle of radius coming back to the 
same point after traversing a circle of radius The fluctuations about this path are taken 
into account by letting 



[0{*M{s)] 



7T 27T S 

- + 5e(s), — + 6<j ) (s) 



(24) 



where 56(s) and 5<f)(s) represent the deviations from the extremum path on the unit sphere 
expressed in terms of angles. Expanding the bending energy of Eq. ( 120]) correct up to second 
order in the fluctuations S9(s) and 5(j)(s) gives 

E bend = (^) 2 + (f + dJ ^) 2 - (f ) 2 Sn*)} ■ (25) 

We expect this expansion to be a valid approximation near the ring configuration, if the 
deviations from the circular configuration is small. 
We expand fluctuations as 

oo 

cos (-j-) (26) 

n=0 
10 



and 

oo 

W(«) = j^W B oos(^). (27) 

n=0 

In terms of these modes, the bending energy of the chain is given by 

oo „ 

E bend = f- L {-^ 2 56l + £ (n 2 - 4)vrW + £ n 2 vr 2 (5^ - + £ nV 3 ^ 2 }. 

n=l n,odd n,even 

(28) 

In the above, the bending energy is independent of the modes 89 2 and <50o- These give two 
of the three rotational degrees of freedom of the ring polymer in space. 8<fio corresponds to 
the rotation about the Z axis, while S62 corresponds to rotation of the ring about X axis. 
A fluctuation of the form <5# 2 s sin(27rs/L) leads to rotation about Y axis. The value of 89 2s 
(the amount of of rotation contained) in an arbitrary 89{s) may be found from 

L 

59 2s = \\ 66(s) sin (^) . (29) 



L J w V L , 


Using Eq. ([27]) in Eq. ([29]), one gets 

69 2a = 2^2a n 5d n , (30) 

n 

where 

a n = r 2~ 4 n , if n IS odd 
n (n 2 -4)7r' r<Q\ 

= 0, if n is even. 

While evaluating the partition function one must avoid integrating over the rotational modes, 
since within our approximation scheme, these modes would cause the partition function to 
diverge. One can remove these rotational degrees of freedom by inserting the product of delta 
functions 5(5<f)o)5(562)5(562 S ) to the functional integral and then taking the contribution of 
the rotational modes to the partition function into account by explicitly putting in the factor 
87r 2 . Then the probability for the loop formation is given by 

G(0,L|0) = §J / D59(s) J D5<f>(s) exp 8(5<j> )8(89 2 )S(89 2s )S(R). (32) 

R is the end to end vector for the polymer chain and Zr is the partition function for the 
polymer, approximated by that appropriate for a semi-flexible rod of length L (see Eq. f[T2l ). 
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B. The asymmetry in the three directions of motion at the ring geometry 

In Ref. we derived the expression for the end-to-end vector R, by expanding the 
components of u(s) as a Taylor series up to first order, which is 

R = L(iJ2 a n 64> n -i 6 -^ + k59 ) . (33) 

\ n odd / 

The components of R are given by, R x = LY^odd^^n-, R y = — f^02 and R z = L59 . 
Thus R x can be changed by varying the value of 5<f) n s for odd n. It can be easily seen from 
Eq. ( BED , that increasing 5<ft n s with n odd decreases the bending energy of the chain towards 
a minimum at 5<p n = 8/n 2 ir (n odd). Using this value for 5(j) n (n odd ) one gets 

OO g 

R x = lY] a n — = L, (34) 

n,odd 

since Yln°odd a n/ n2 — tt/8. Hence, this value of 5(p n corresponds to a rod lying along the X 
axis. Therefore, the ring is unstable along R x and the bending energy along this direction has 
the minimum at R x = L. Also, increasing 58q decreases the bending energy and therefore, 
R z is also unstable. This is because when R z is increased the ring changes into a helix, which 
has less curvature and therefore less bending energy. (Note that we do not take torsional 
energies into account in this analysis). But unlike R x = 0, R z = is a maximum. It should 
be noted that R z = L corresponds to the rod and should be a minimum, but our analysis 
does not reproduce this. So, the instability along R z is only near the ring, where our analysis 
is valid. Unlike R x and R z , R y is stable, since the bending energy of the ring increases when 
502 is increased. Thus, the motion in R x , R y and R z directions are energetically different. 



C. Bending energy of the rod 

For a semiflexible chain, the minimum energy configuration is the rod. On the unit sphere 
representing the tangents this means that the random walker stays at the starting point. We 
take this point to be {0(s), 4>{s)) = (vr/2, 0), which corresponds to the rod lying along the X 
axis. Unlike the great circle, the straight rod is an exact minimum energy configuration. In 
this case, the fluctuations can be incorporated by letting 

(35) 
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[6(8),<K8)] 



TV 



+ 59(s),5(j>(s) 



The bending energy of the rod correct up to the second order in fluctuations is then given 
by 

o ^ ' 

Using the expansions Eq. (126!) and Eq. (127|) . one gets 

2JV 
n=0 

Unlike the ring, the rod has only two rotational degrees of freedom. From Eq. ( 1371) . it follows 
that these are the modes 5(fto and 59q. For the convenience of bookkeeping, we assume that 
the number of 8(p n and 59 n modes and both equal to 2Af, with n = 0, 1, 2...{2Af — 1). Thus 
there are AM modes in total, with M — > oo. 



IV. THE NORMAL COORDINATES AND THE GREEN'S FUNCTION 

In this section, we use the approximation scheme for the bending energy described in 
Sec. IIH Al to analyze the dynamics of loop formation. The approximation of Eq. ( 1281) 
is valid only near the most important path corresponding to the loop, since the fluctua- 
tions about this path are assumed to be small. The time evolution of the chain may be 
described by the multidimensional Green's function Gq(^ ,t\^ Q ), where = (<&', &') with 
$t = (5(f) ,5(f)2, ■■S(f>2M-2, S<f>i, 803, ■■8 ( t ) 2Af-i) (note that we have separated out the even and 
odd modes) and 0^ = (S9 Q , 59 2 , ••<56 , 2Af-2, SOi, S9 3 , ■ ■S9 2 j^-i)- The superscript f stands for 
transpose. The bending energy of the polymer near the ring configuration is given by Eq. 
( 1281) and therefore, Go for configurations close to the ring may be obtained by solving the 
corresponding equations of motion. Because of our approximation for the energy, the Green's 
function so obtained is not valid for large R. Yet, the sink-sink correlation function, T>(t) 
of Eq. (fl2l) . may still be evaluated, since the sink function <S(R) is nonzero only for very 
small values of R. This of course, is approximate. The function Go may be found by solving 
the equations of motion of the polymer near the loop configuration. The angle coordinates, 
^ are not normal coordinates, since the kinetic energy of the chain has terms that couple 
these (see Appendix [A]) . As a result, the equations of motion of the chain in terms of them 
are coupled. This coupling may be avoided by working with the normal modes, which may 
be found by solving the corresponding eigenvalue problem (Eq. fpf6]) ). Then the dynamics 
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of the chain can be reduced to the dynamics of a particle in a multidimensional harmonic 
potential. The Green's function Go is obtained as a product of the one-dimensional Green's 
functions corresponding to each of the normal modes. 



The Hamiltonian and the normal modes 



The kinetic energy of the polymer in the center of mass frame is 



T= P - 



L 

f 


dr(s) 




dt 



dt, 



(38) 



where r(s) is given by Eq. ( 122|) . Near the rod ring configuration, the kinetic energy of the 
polymer in terms of the Fourier modes 59 n and 5<ft n may be written as 



R 



(39) 



The dot in represents differentiation with respect to time. The subscripts R{L) in 
^rI^l) are used to indicate that these are deviations measured from values appropriate 
for the rod (loop) geometry. Tr is the kinetic energy matrix appropriate near the rod 
configuration. It has a block diagonal structure, having no matrices connecting the 9 and 
modes. Even within the 9 and modes, odd and even modes are decoupled. Hence Tr 
may be written as 



T 



R 



rj-i0e 

R 














L R 














*-R 

















Detailed structures of the T matrices is given in Appendix In a similar fashion, near the 
loop configuration the kinetic energy is given by 



Tr 



(40) 



The angles ^ L 



($i, f ) with & L 



(50o, 502, ~8<pwr-2, <50'i, 503, ••502 A r_i), where 



50 n — 8/(n 2 7r). Like Tr, Tl too has a block diagonal structure, with the blocks given by 
the matrices T^ e , T^°, T^ e and T^°. The forms of these too are given in Appendix |A] Note 
that these matrices have no length (L) dependence. It is found that (see Appendix [SJ the 
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modes of odd and even n decouple. One may rewrite Eq. ff28l) as 



^bend 



IL 



-M 2 56 2 + J2 {n 2 - A)7i 2 Se n ' 



n even 



K 

IL 



n,odd n,even 



+ 

(41) 



or as 



£ 6e „d = ^*[V L * L . (42) 
Thus the total energy of the polymer molecule near the loop configuration is 

E = ^[r L 9 L + ^*[V L * L . (43) 

Like Tl, Vl too are block diagonal. The matrices of which Vl is composed of are 
V^ e , V^°, V e L e and V^°. Each one of them is diagonal and have matrix elements given 
by (V^°) nm = 5 mn n 2 n 2 /2 = (V^ 6 ) nm ; (V'j')oo = — 47r 2 and all other matrix elements being 
given by (V L ) nrn = S mn (n 2 - 4)vr 2 /2. 

We will choose the sink function S(r) as a function only of the end-to-end vector R 
(see next section). The dynamics of the closing process must be unaffected by the spatial 
rotations of the polymer. Hence, the sink-sink correlation function, T>(t) of Eq. ffl2l is 
independent of the rotational modes 5<j)o,582 and 682s- From Eq. ( J33l) it follows that S6 n 
modes with odd n do not contribute to the end-to-end separation R of the polymer. As R 
has no dependence on the odd 58 n modes, they are irrelevant for the dynamics of closing 
process. 

We define Y by 

T 1/2 * L =WY, (44) 
where U is to be defined below. Then the energy becomes 

E = ^yVzYY + ^-YWT- 1/2 V L T- l ' 2 UY. (45) 

Taking hi to be a unitary matrix, which diagonalizes T~ 1 I 2 VlT _1 / 2 to give the diagonal 



matrix /C, as 



WT 1/2 V L T- 1/2 U = K, (46) 

we get 

2 2L 
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E = ^Vy + ^rYtKY (47) 



with 































K ee 

















(48) 



The block diagonal structures of T and V imply that U also has a block diagonal structure, 
with matrices U^ e ,U^°,U ee and U e ° occuring along the diagonal. The energy may be 
written in terms of the components of Y and K as 

3 4Af AM 
n=l n=l 

Note that k n are not dependent on the length L of the chain. Of the modes Y, Y n , with 
n = 3Af + 1 to 47V arise from odd 58 n . The end to end vector R has no dependence on 
them. Hence these Y n play no role in the dynamics of loop formation, occuring near the loop 
geometry. Therefore, we focus on the remaining modes. We write the remaining normal 
coordinates Y/=(Y"i, F 2 , Yzm) as 

Y 7 = (x,y,z). (50) 

i.e, Y i = (xi, X2, --xm, yi, V2, ■ ■■■VM, zi, z%, ....zjs). Within x (y or z), we take the modes to 
be arranged in the order of increasing eigenvalues, and we label them as k xn (k yn or k zn ), 
with n varying from 1 to N '. x n are the normal modes corresponding to the modes the S(p' n 
and these are all stable modes, as may be inferred by looking at the expression for energy of 
Eq. (|28|) . y„ corresponds to the even S(j) n s. S(f> is a rotational mode and correspondingly, 
k y i = 0. Of the even 58 n modes, one is unstable, viz., the one that corresponds to 58o- It 
leads to separation between the two ends in the Z-direction and is unstable, as we already 
discussed. 682 is a rotational mode and would give us a zero eigenvalue. Thus we have k z i 
negative and k Z 2 equal to zero. Note that the eigenvalues k n have no dependence on k or L 
and hence are universal numbers. Their values up to n = 10 are given in Table [B 

The end-to-end distance R may be expressed in terms of the normal coordinates Y/. The 
rc-component 

R x = L }] a n S4>n = Ly j a n 8(j)' n H X^HF' ( 51 ) 

n n n 

which on using a n = — 4/(n 2 — 4)n {n odd) becomes 

R x = a nH'n + L - ( 52 ) 

n odd 

16 



Now, in terms of the normal coordinates (i.e., inverting Eq. 



^2 (t^°) ^ 2 



m=l 



x m . Hence, 



r, x — l — l f n 



71=1 



with 



Similarly, 



where 



fn — 02fc_l 
k=l 



fcn 



AT 



n=l 



fjn 



X ; J 2n 



(53) 



(54) 



(55) 



(56) 



Note that y\ is a rotational mode and the corresponding g% would be zero. Therefore, the 
above sum may be modified to 

AT 

(57) 



n^2 



The z component 



with 



R z = Lj2hr 



n=l 



h, 



In 



Again, Z2 being the rotational mode, this may be written as 

Af 

R z = L ^ h n z n . 

n^2 



(58) 



(59) 



(60) 



B. The equations of motion and the Green's function 



1. The equations of motion 



The equation of motion of the polymer in a dissipative environment is given by 

5E[v{s)\ 



pr + pjr + 



St{s) 



CM), 



(61) 
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n 


k 




h 


g n /ky n 


h 




1 


298.54 


3.7839£ - 


02 


0.0 


0.0 




-1241.1 


-2.5996£ - 


- 02 


2 


7.1266£ + 03 


1.5524B - 


04 


1.4131^ + 03 


1.1665£ - 


02 


0.0 


0.0 




3 


6.3197£ + 04 


1A907E - 


06 


2.4549£ + 04 


7.9245£ - 


04 


2.8705£ + 04 


5.0074S - 


04 


4 


2.5382£ + 05 


6.3190£ - 


08 


1.3492S + 05 


1.3331E - 


04 


1.5259S + 05 


1.0042E - 


04 


5 


7.0125S + 05 


5.8587£ - 


09 


4.3631£ + 05 


3.9529£ - 


05 


4.741LE + 05 


3.2983£ - 


05 


6 


1.5655£ + 06 


8.6906£ - 


10 


1.0696S + 06 


1.5805£ - 


05 


1.1339^ + 06 


1.3963£ - 


05 


7 


3.0455£ + 06 


1.765&E - 


10 


2.2137£ + 06 


7.5909£ - 


06 


2.3108£ + 06 


6.941LE - 


06 


8 


5.3826£ + 06 


4.5U2E - 


11 


4.0872£ + 06 


4.U20E - 


06 


4.2234£ + 06 


3.8728£ - 


06 


9 


8.8695£ + 06 


1.389LE - 


11 


6.9516£ + 06 


2.4991£ - 


06 


7.1330^ + 06 


2.3728S - 


06 


10 


1.3893S + 07 


5.1570£ - 


12 


1.1123^ + 07 


1.6576£ - 


06 


1.1355^ + 07 


1.5914E - 


06 



TABLE I: The dimensionless ei 



genvalues k n and the corresponding values of fn/k xn , g n /k yn 



and 



where J5[r(s)] is the energy functional of the chain and ({s,t) is the stochastic force act- 
ing on the s th segment of the chain. (,(s,t) is assumed to obey < ((s,t) >= and 
< £(s, t)C(s', f) >= 2kBTp r y5(t — t')5(s — s'). In the over-damped limit, one may write 
Eq. flEQ as 



P7r + 



5E 
Ml) 



C(s,t). 



(62) 



Through the use of a system-plus-reservoir model, this equation can be equivalently ex- 
pressed in the angle coordinates as (see Ref. [l|) 



(63) 



and 



(64) 



Equations (1631) and ( 1o4l) represent sets of coupled first order differential equations. For the 
ring, we can express them in terms of the normal modes. In terms of the normal modes Y n , 
equations (163]) and (}6"lj) represent a set of independent one dimensional Langevin equations, 



Yn + . Y n £ n (t) , 



(65) 
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where ( n (t), is a white Gaussian noise with < ( n (t) >= and < ( n (t)( m (t') >= 
/(pL 3 )5(t — t')5 mn . From Eq. (1651) . it follows that the relaxation time of each mode 

~ L 4 . 



2. The Green's function 

Eq. ( 1651) describes a particle of mass pL 3 , subject to friction 7 in a one-dimensional 



harmonic potential, Kk n Y^ /2L. The Green's function for it is given by 35l. 1361] 

-1/2 

(66) 



G n (Y n X;t) = {^[l-exp(-2t/r n )] 



cxp 



Kk n (Y n -C X p(-t/Tn)Y^) 2 

' 2XA s T[l-exp(-2t/T n )] 



with r n = To/k n with tq = pL 4/ -f/n. G n is the conditional probability to find the particle 
at Y n at time t, given that it was at at t — 0. Eq. (1661) is valid for both positive and 
negative k n |36j. 



V. THE EQUILIBRIUM DISTRIBUTION 

In this section, we derive an approximate equilibrium distribution function for the semi- 
flexible polymer near the loop configuration in terms of the angle coordinates. The partition 
function of the polymer is 

Z = J dr J dpexp(-/3#[r,p]). (67) 

In terms of the angle coordinates, the partition function is given by 

Z = J d<S> J rfp*exp(-/3if[*,p*]), (68) 

where p^ are the momenta conjugate to the angle coordinates ^ . In the integral, the 
configurations that contribute the most are the ones near the rod configuration. For such 
configurations, the energy is given by (see Eq. (1371) ) 

E = ^*Wr*R + (69) 
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where Vr too is block diagonal, consisting of 



V 



R 



V 



<j>e 
R 







o v£° 





























(70) 



The matrix elements of each of the matrices on the right hand are given by (V^ e j = 

V / ran 

= {K e ) mn and (V$) = ^fd mn = (V£) . The momenta conjugate to * R is 

V / ran 

P* H = pL 3 Tft 1 '*S? R and hence the Hamiltonian is given by 



(71) 



The partition function of the rod can be evaluated now, using equatioins (J55j) and (I7T1) . 
polymer given in Eq. ( 1681) near the rod can be evaluated now. The rod has two rotational 
modes, which are S4>q and 58 . Integrating over them would give a factor of Air. Performing 
the integration over the remaining 5<fr n and 58 n and integrating over all the momenta, give 



Z R = 4n(2npL 3 k B T) 2M (detT R ) 1/2 {2nk B TL/K) 2Af -\det'V R )- 1/2 . 
1 



^L^-^^K^p 2N 



■ R\ 



R 



(72) 
(73) 



((2^-l)!) 2 ' 

The prime on the determinant in det'Vi? indicates that the zero eigenvalues (rotational 
modes) are excluded. Note that we use |T^| to denote |T^, 6 ||T^°| and \T e R \ to denote 
|T^ e ||T^|. For configurations close to the loop, the Hamiltonian can be approximated by 



(74) 



This can be used to calculate the equilibrium distribution near the loop conformation. In 
particular, we are interested in the probability of contact between the two ends at equilibrium 
G(0,L) defined by 



G(0,L) 



1 

Z 



dp* / exp (-/?#) 5 (R). 



(75) 



5(R) in the above ensures that the two ends of the chain are in contact. Our strategy in 
the calculation is as follows. The major contribution to the partition function comes from 
rod like conformations. Hence we approximate Z = Zr. Near the loop geometry, we can 
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use the approximation H = Hl and perform the integrals over the angles, with rotational 
degrees of freedom easily accounted for. Thus 

G(0, L) = J dp* L J d* L exp(-f3H L )5(R). (76) 

The integrals over momenta are easy to perform. There are three rotational degrees of 
freedom, which can be removed by inserting 5(5(f>o)5(592c)S(S62 S ) into the integrand, and 
their contribution accounted by introducing a multiplicative factor of 8n 2 . Then the 
integrations can be performed, one by one, after using the integral representation for 5(R) 
and using Eq. (|33|) for R. Thus 

G(0,L) = ^~f dp* L J d^ L exp(-i3H L )6(R)5(S(j )o )5(5e 2c )8(Se 2s ). (77) 

The details the calculation are given in the Appendix [B] and the result is 



G(0,L) = j=J (78) 



L 



3L 5 



with 



rp<p 

L 



det(T^). Putting in numerical values (see Eq. (IBldO ). we find 

G(0,L) = 1522.06e"^^--^ (79) 

Lj 



a form that is in agreement with the results of Shimada and Yamkawa [26[|. Note that the 
persistence length of the chain, l p — (3k. We give a compartive plot of the our function and 

their function G SY (0,L) (= 896.32e — +0 - 246L / l ?(l p /L) 5 ) in Fig. El It is clear that there is 

fair agreement between the two. The value of L at which the maximum occurs is L = 2.63/ p 
in our G(0,L), while it occurs at L — 3.37/ p for the results of Shimada and Yamakawa [26 ]. 

It is interesting to ask how the L~ 5 term in Eq. (!79|) comes about. The Dirac delta 
function 5(R) contributes L~ 3 . G(0, L) would have the ratio of the partition function 
for the loop conformation to that of the rod, and this contributes a factor of {(3k/ L) 1 ^ 2 . 
Further, the fact that the potential energy depends on k/L term causes three factors of 

1 /2 

{(5k/ L) (from the three components of R contained in the probability density at i?j = 0, 
with % = x, y or z and comes from the fact that the larger the value of L, the broader the 
distribution of the Rj). These multiply together give the factor (3 2 k 2 /L 5 in Eq. f!79|) . 
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0.10 



E 

Q 0.05 



0.00 



- Our result 
SY 




2 3 4 5 6 7 
L (Dimensionless) 



FIG. 2: Comparison of our result for G(0, L) (full line) with the result of Shimada and 



Yamakawa 



26l ] (dashed line). 



VI. THE TIME FOR LOOP FORMATION 

We now evaluate the average loop formation time using the approach outlined in Section 
HTl The quantity v eq is 

v eq = (S(B)) (80) 

with (3(B)) defined by 

(S(B))) = ^Jd* J dp 9 exp(-PH[i&,p 9 ])S(R). 
Following our discussion in previous Section, we approximate it as 



(S(B)) = ±- J d* L J dp* L exp(-/?tf L [*,p* L ])S(R) 
It can be easily evaluated, following the methods of Appendix [B] The result is 



(S(R)> = 16 



Ly/L 3 - 4vr 2 /?r/ 2 K V /L 3 + 8vr 2 /?r/ 2 K^3L 3 + ^ 2 (3r] 2 K 



(81) 



(82) 



(83) 
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A. The sink-sink correlation function 



The essential step in finding the average time of loop formation is to calculate the sink- 
sink correlation function, Eq. (fT21) . The sink-sink correlation function can be written in 
terms of 

V(t) = f d* f d*'S(R')G (*',*;t)S(R)P eq (*). (84) 

Go(ty', t) is the propagator expressed in terms of \I/ and obeys the condition 
G (*', *; t)S{R) -> - *)) as t -»• 0. The P eq {^) in the above is given by 

i% W = ^Jd* Jdp* exp(-/3# [*, p*]). (85) 

In the spirit of our previous discussions, we approximate P eq {^) near the loop configuration 

as 

Pe q {*) = Y~ R \ d P* L exp(-/3# L [*,p* L ]). 
The sink-sink correlation function may be written as 

V(t) = (S(R)}C(t) } (86) 

where 

C(t) = y d» y ^ / l S(R')G (* / , *; t)5(R)P eg (*). (87) 

C(t) can now be approximated as 



/ d* L J dV L S(R')G (V L ,* L ;t)S(R)exp (-(3^[v L ^ L 

C(t) S C Q (t) = 7 ■ ^ '-. (88) 

Jd* L <S(R) exp (-/3^*t.V L * L 



Note that we use the subscript "a" to denote the approximate value of C(t). The above 
integral may be re-expressed in terms of the normal modes Y as 

/ dY' J dYS(R')G (Y', Y; gg(R) exp (-g^ytgy) 
^ j / c/Y<S(R) exp (-/?^Yt/CY) 

(jo(Y', Y; £) is the propagator expressed in terms of the normal modes Y. The Jacobians 
associated with the transformation in the numerator and denominator cancel out (Note that 
/ dY'G (Y\ Y; t) = 1). With the above form the sink function, C a (i) can be evaluated to 
obtain (see appendix ICl) 

C a (t)=C x (t)C y (t)C z (t) (89) 
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with 



C x {t) 



2 K 



2iry/{S x {0)L 3 + W") 2 ~ L 6 S x (ty 



and 

On using these, 

C a {t) = 



Cy(t) 



C z {t) 



2 x/2tt 3 / 2 y/{S y {0)L 3 + (3ri 2 K) 2 - L e S y {t) 2 

VP^S z (0)L 3 + (3 v 2 k 
V^tt 5 / 2 ^{S z {0)L 3 + /V«0 2 - L^S z {t) 2 ' 



(90) 
(91) 

(92) 



\ 



3 






2 


rp<j> 

L 





\L 2 /3k 



Sx(O)L'^+0 v 2 k (S x (0) + S x (t))L^+l3 v 2 K + SL^+Stt' 2 ^,, 



128 



(93) 



\/(L 3 - 2ttW«) (^ 3 + 8ttWk) (3L 3 + 8ir 2 Pri 2 K)y/S x {0)L 3 + Pri 2 Ky/S y (0)L 3 + (3r] 2 n^/2S z (0)L 3 + 2/?r/ 2 , 

T^Op + /fy 2 «) 2 - L 6 ^(t) V(S y (0)L 3 + 2 - £ 6 ^(t) V(^(0)L 3 + /V«) 2 - Z,6 5z ( t )T 

In the above 



n=l rc=2 717^2 



A/' 



(94) 



with t = t/r . 5^(0), S y (0) and 5 Z (0) can be evaluated exactly (see Appendix ID]) . Using 
their values given by equations (1D2I) . (1D3h and (1D6h . Defining 



/(f) = S x (t)/S x (0),g(t) = S y (t)/S y (0);a,nd h(t) = S z (t)/S z (0) 



(95) 



we get 
C a (t) -- 



2(3k 



71 



3/2 i 



+8tt 2 /3t7^ k 3/(t)i a +3i a +8 tt^/Stj 



(96) 



v/L 3 - 4 n 2 f3ri 2 Ky/L 3 + 8n 2 (3r] 2 n^3L 3 + 8vr 2 /?r/ 2 K 



(3L 3 + 87r 2 (3r] 2 Kj z - 9L e f(t) 2 J (L 3 + 8vrW/t) z - L 6 #(f) 2 J L 6 h(t) 2 - (L 3 - 4ttW^ 2 

The values of f n /k xn , g n /k 2 n and h n /k 2 n are given in Tabled! Use of equations (l8~3l) and 
( 1961) in Eq. (1861) leads to an approximation for T>(t) which we denote by T> a (t). 



V a (t) = (S(R))C a (t). 



(97) 



Since the eigenvalue k 2z is negative (see Table [I]), the term g(t) has a term that diverges 
exponentially as t — ► 00, making C 2 (oo) = 0. Hence V a (oo) is zero. This is due to the 
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instability along R z , which causes the correlation function to vanish at long times. Hence, 
T> a (t) given by Eq. fl97|) is a good approximation to T>(t) at short times. For long times 
V(t) will approach v 2 eq . Hence V(t) = V a (t) is a valid approximation only at short times. 
For for long times, T>(t) should be equal to v 2 . Hence, it follows that the actual correlation 
function may be approximated as T>(t) = T> a (t) + v 2 . Using this in Eq. (ITTj) one gets 



T = 4" / T> a {t)dt 



(98) 



We note that /3k = l p is the persistence length of the polymer and that 7 has dimensions 
of time, and use these as units for length and time. Then, the expression for r becomes 



\ 



3 


It 




64 




7T 9 



L g / 2 W(L) 



(99) 



with 



W(L) 



-1/2 



f(ty 



[1 + 87r 2 r7 2 /3L 3 



-1/2 



(100) 



1 - 



h(tf 



-1/2 



1 + 8ttV/^ 3 ) / V(1-4ttV/^ 3 ) 
The functions f(t), g(t) and h(t), on evaluation are found to be given by 



fit) = 0.995873e 



-298.541* 



+ 0.00408587e 



-7126. 56t 



(101) 



g(t) = 0.921041e- 1413 ' 12 ' + 0.0625692e" 24548 - 8 * + 0.0 1 05 2 56 e - 134918t 

+0.00312111e~ 436308 i + 0.00124791e- L06959xl ° 6 *. (102) 

hit) = 1.02629e 124L14 * - 0.0197683c- 28705 ' 1 * - 0.00396455 e - 152588 - i 

-0.00130213e- 474114 -*. (103) 

Terms that make no significant contribution to f{t), g(t) and h(t), as their exponents are 
large and the coefficients small have been neglected in the above. Using Eq. ( |99l) . we have 
calculated the average time of loop formation r as a function of the length L, and the results 
are given in figures [3] and HI The results are dependent on the sink width rj. The full lines 
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in these figures give exact values of r for various values of sink size 77. The values of r for 
each value of 77 were fitted with an equation of the form A exp(E a / L)L n and the values of E a 
and n, as well as the value of the length at which r is a minimum (L m ) is given in Table HT1 
The curves that result from the fitting are shown as dotted lines in the figures. It is seen 
that the functional form r = Aexp(E a / L)L n reproduces the data well, with the exponent 
in the prefactor n being dependent on 77 and varying from 4.5 to 6. For small 77, n is close 
to 6. Further, the loop formation time r becomes longer and longer as 77 is made smaller. 
In fact, as may be seen by from Eq. (I99j) . the time diverges like I/77 as 77 — > (see Eq. 
(11051) ). This is not surprising, because as 77 — > 0, the diffusional search for loop formation is 
for a smaller and smaller volume in space. In fact, on looking at Eq. (I99p and remembering 
that 77 actually serves the purpose of a small length cut off, one would have expected the 
dependence to be approximately 

r = Aexp(E a /L)L 9 / 2 . (104) 

However, there are two reasons that lead to the observed dependence on 77. (1) The expo- 
nential term in Eq. fllOOj) is dependent on t and makes the L dependence change from the 
simple form of Eq. (11041) and (2) the L dependence of the terms inside the square roots 
in Eq. (1 1001) . From the functional forms of /(£), g(t) and h(t), it is clear that due to the 
presence of h(t), the integrand decreases rather rapidly. For t ~ 1/ot,, where 07, = 1241.4 
(see Eq. [TOT]) ( {1 _^2 /L sf ~ l) - H*)' 1 = (1/1.02629)6^* . On this time scale 

(t ~ l/w h ), one can approximate (l - j^^fjj^) ' = 1 and (l - j^^fj^jf) 
= (3L 3 /167r 2 ry 2 ) 1,/2 . Hence for small values of rj, the integral may be evaluated approximately 
to get 



7~close 



3 








3207,77 \ 






7T 11 



L 6 e^r. (105) 



Thus for small 77 the loop formation time behaves like ~ L 6 e sl as seen in the Ta- 

1—1 / f(t) 2 \~ 1/2 

ble |HJ For not so small values of 77, one expects ( 1 — 2 ' 2 3 a ) — 1 and 

\ (l+o7r rj I AL ) J 

(l - / -1 ft IT - v A = 1. This leads to n ~ 9/2. 

Physically, the above results are easy to understand. The rate of the loop formation may 
be written as ~ V(R X < 77, R y < 77, R z < r/)x frequency factor xZ L /Z R . In this, V(R) 
is the probablity distribution function for the end-to-end vector. For a semi-flexible chain, 
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Exact 




1 1 1 1 1 1 1 1 1 1 1 1 1 r 

1 2 3 4 5 6 7 



L (Dimensionless) 

FIG. 3: The time of loop formation, r as a function of L, for different values of the width r/. Units 
are chosen such that both r and L dimensionless. The full curves are the computed results. They 
are well represented by the functional form AL n exp(E a /L) as may be seen from the figure, where 
we have represented them by dotted lines. The parameters that result from fitting are given in 
Table IE 

the frequency factor ~ L~ 4 . In the limit i] — > 0, V(R X < 77, R y < 77, R z < 77) ~ {(5k/ X) 3 / 2 as 
seen earlier and Zl/Zr ~ ((3k/L) x I 2 . Therefore, the preexponential factor of the rate has 
((3k) 2 /L 6 dependence. On the other hand, as one increases the value of 77, for sufficiently 
large r], V(R X < 77, R y < 77, R z < if) = 1, leading to rate of the form ((3k) 1 / 2 / L 9 / 2 (37)]. 

B. Numerical Results 

We now consider loop formation of double stranded DNA, which has a persistence length 
of 50 nm. The sink function, defined by Eq. f|T3|) and f|T4l) has a width equal to 77 which may 
be taken as 2-3 A. Then the dimensionless rj would have a value of roughly 1/200 and 
this would correspond to the lowest curve in Fig. HI with n = 5.9. On the other hand for a 
more flexible chain, with persistence length equal to 2.5 nm, and with a value of r] equal to 
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FIG. 4: The closing time r as a function of L, for different values of the width rj. Units are 
adopted such that all these are dimensionless. The full curves are the computed results. They are 
well represented by the functional form AL n exp(E a / L) as may be seen from the figure, where we 
have represented them by dotted lines. The parameters that result from fitting are given in Table 

m 

1 A, one would have dimensionless rj = 1/25, and this would correspond to the lowest curve 
in Fig. [3] with n = 4.9. The value of L m at which the minmum time is required for loop 
formation, does not depend strongly on the value of rj. Thus it is found to be in the range 
2.2 to 2.4 times the persistence length of the chain (see Table HT|) . 

The dynamics of loop formation in semiflexible polymers was analyzed by Dua and Cher- 



ayil, who found r ~ L n , with n ~ 2.2 — 2.4, with n approaching 2 in the flexib 
is obviously valid in the longer chain limit. On the other hand, Jun et al % 



e limit. This 
have studied 



the region where the length of the chain is a few times the persistence length. The assumed 
the two ends of the chain to execute random walk with a constant diffusion coefficient, and 
found that there is a length (L m ) at which r is a minimum. Their analysis used accurate 
results for G(0,L) and lead to somewhat larger value for L m (3 — 4). On the other hand, 
we have studied the dynamics in detail, using a multimode approach. We get expressions 
for G(0, L) and r which have extrema at lower values of L, this being a result of the use of 
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V 


A x 10 6 


E a 


n 




0.04 


2.800 


10.247 


4.920 


2.19 


0.02 


1.823 


11.982 


5.261 


2.32 


0.01 


1.627 


12.890 


5.51 


2.37 


0.005 


1.740 


13.421 


5.74 


2.37 


0.0025 


2.268 


13.64 


5.90 


2.33 


0.00125 


3.670 


13.65 


5.99 


2.28 



TABLE II: Fitted Parameters: Values of A, E a and n and L m (the dimensionless length at which 
r is a minimum). 

approximate expression for the bending energy. 



VII. SUMMARY AND CONCLUSIONS 



In this work, we have presented a detailed mult i- dimensional analysis of the loop for- 
mation dynamics of semiflexible chains. The reverse process, the opening of the loop was 
studied in a previous work [1] , where we developed an approximate model for a semiflexible 
chain in the rod limit. In this model, the conformations of the polymer are mapped onto 
the paths of a random walker on the surface of a unit sphere. The bending energy of the 
chain was expanded about a minimum energy path. This model was shown to be a good 
approximation for the polymer in the rod limit and provided an easy way to describe the 
dynamics. Use of this model led to opening rates of a semiflexible polymer loop formed by 
a weak bond between the ends, in Ref. lj, we calculated the opening rates for a Morse 
type interaction between the ends of the polymer as a function of the contour length of the 
chain. In this paper, we analyzed the loop formation dynamics using this model and thus, 
presented a rather complete theory of dynamics of formation of semiflexible polymer loops. 

The dynamics was described using the formalism by Wilemski and Fixman, which de- 
scribe the intra-chain reactions of polymers as a diffusion-controlled reaction. In this formal- 
ism, the reaction process is described using a sink function. For an arbitrary sink function, 
exact results are not available and hence, WF introduced an approximation called "closure 
approximation" . In this procedure, the closing time can be expressed in terms of a sink-sink 
correlation function. To calculate this sink-sink correlation function and thereby the closing 
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time, one needs to know the Green's function of the chain and the equilibrium distribution. 
We calculated the Green's function of the chain through a normal mode analysis near the 
loop geometry. This normal mode analysis could be performed independently of the rigidity 
(«) and contour length (L) of the polymer, leading to a set of eigenvalues that are universal. 
An approximate equilibrium distribution for the polymer near the ring configuration was 
given. As the sink function vanishes for large values of the end-to-end distance R, sink- 
sink correlation function has contributions mostly from the dynamics of the polymer near 
the ring configurations. We calculated this approximate sink-sink correlation function for a 
Gaussian sink through a transformation of variables into normal coordinates. 

We then obtained loop formation time (in dimensionless units), r for different contour 
lengths of the chain. We found that r ~ L 9 / 2 W(L), where W(L) is an integral that could be 
performed numerically. Numerical calculations lead to the result that r = AL n exp(E a /L), 
with n varying between 9/2 and 6. r was found to have a minimum at L m j n — 2.2 to 2.4 
which is to be compared with the values 3.4 obtained by Jun et. al. [3] by a simple one 



dimensional analysis and the value 2.85 of Chen et. al. 



29J found through simulations [29] 
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APPENDIX A: THE KINETIC ENERGY OF THE RING AND THE ROD 

The kinetic energy Eq. (1551) of the polymer may be evaluated using Eq. (T2"2"l) . We take 
r cm = 0, so that the ring is described in the center of mass frame and the translational 
degrees of freedom are eliminated. 

1. Matrix elements for the Loop 

The kinetic energy matrix elements of the loop are given below. 
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a. Odd 4> Modes 



For odd 8<f) n modes one has 

(T^ )rrm t n 5 mn 16t n t m , (A-l) 



where 

4 + n 2 

tn = 2(-4 + n 2 )% 2 - (A2) 

b. Even <p Modes 
For even 5(f) n , one has 

(T^ ) mn = t n 5 mn (A3) 

with n, m ^ 2. For m,n = 2 one gets 

flt% - =?j±£ (A4) 

and 

^L)2n = - 2( _ 4 + n2)27r2 = (Tf)n2 (A5) 

c. Odd 9 Modes 

The kinetic energy matrix corresponding to the even 9 modes have the following form. 
For n,m/0 

(T^ )nm Qn$mn (A-6) 

where 

9n = 7zkl5- ( A7 ) 



2n 2 7T 



For the zeroth mode 



(Ti°)oo = ^ (A8) 



and 



(T£V = ~J- 2 = (T 9 L ) n0 . (A9) 
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d. Even 9 Modes 

The kinetic energy matrix corresponding to the odd 9 modes have the following form 

C^-L )nm = Qn^mn ~ d n d m , (AlO) 

where 

(All) 



w 2n 2 7r 2 
and 



d n = (A12) 



n 2 Tc 2 



2. For the Rod 

In the case of a rod, the matrix elements are identical for <p and 9 modes. They are given 

by 

(Tr )mn = (T^°) m „ = t n 5 mn — 16t n t m , (A13) 

where 

'« = ss?- < A14 > 

For modes with even n one has 

(T'ft ) mn — (T^) mn = t n 5 mn (A15) 
with n, m ^ 0. In this case, the zeroth mode is couples to the other modes. Thus one gets 

(T?)oo = (T 9 R e )oo = \~ 2 (A16) 

and 

(If V = (TtV = J~ 2 = (Tf)no- (A17) 
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APPENDIX B: THE EVALUATION OF G(0,L) 



We now give details of the evaluation of G(0,L) . We perform the integral in Eq. (|77|) 
and substitute the value of Z R from Eq. (1721) to get 



G(0,L) 



8vr 2 



d v * L / dy L exp(-pH L )5(R)5(5<f> )5(59 2c )5(56 2s ). 



7T2 



3-4Af 



2AT-1 



\ 



(Bl) 
(B2) 



where G^°G^ e G e °G ee are defined and calculated in the following. is the contribution 

from the odd <ft modes to G(0,L) and is defined by 



G 



° = n / d<W « 

n,odd 



exp 



4L 



(B3) 



We put 5(56 2s ) = (1/2tt) / dp exp(ip56 2s ) = (1/2tt) / dpexp(2pi £ aJ6 n )- With this, 
integrals over with n = 3,5, .. are evaluated and then the one over p, after which one 



can easily evaluate the integral over 59%. The result is 



G e 



8o _ ( ™P 

L 



(l-AO/2 



^T(Af- l/2)r(jV + 3/2) 



(B4) 



T is the Gamma function. G 06 is the contribution from the even # modes to G(0,L) and is 
defined by 



^ = II / * ex P 



AL 



n 2 - A)59 2 n 6(56 2 )6(L6e ) 



The 5(L59 ) comes as the 5(R Z ) part of 5(H) in Eq. (1B1I) . The integrals are easy and the 
result is 

G^ e is the contribution from the even <p modes to G(0,L) and is defined by 



= II / d6 ^exp 



n,even 



f3nn 2 
AL 



n 2 5<p 2 n 5(5^)5^5^/2) 



(B6) 



The 5(5(po) comes as a result of removing the rotational mode 5(po- ( IB II) . The integrals are 
easy and the result is 
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is the contribution from the odd <j) modes to G(0,L) and is defined by 

(3kti 2 



(;■"■ = ] [ 

n,odd 

On evaluation, we get 



G*° 



cxp 



AL 



-n 2 (^) 2 )5(L^a„(^; + 

n,odd 



n 2 7T 



7TK/3 



(l-jV)/2 



exp( 



4tt 2 k/? n 



(B8) 



vf^IV L / T(A/+l/2) " v 3L ' 

The equations ( 1B4I) . ( 1B5I) . ( 1B8[) and ( 1B7[) above combined together with ( 1B1I) and with 
Af — > oo taken, gives 



G(0,L) 



lev^Tr 3 /? 2 



3L 5 



\ 



cxp 



4tt 2 /?k 
3L 



(B9) 



The ratio 



can be evaluated using MATHEMATICA, taking each to be 1000 x 1000 



matrices. It evaluates to 



\ 



R 



6.5083. 



(B10) 



APPENDIX C: THE EVALUATION OF <5(R)) AND C(t). 

The evaluation of (S(R)) and C(t) are similar to the evaluation of G(0,L) carried out 
in Appendix IB! The result for (5(R)) is 



(S(R)> 



\ 



R 



16e 3i3+i^W^7r 3 /5 2 /« 2 



(CI) 



C(t) may be written as a product of three terms, as already seen in Eq. (1891) . We give 
expressions for them the following. 



C x {t) 



flJdY'J dY n G(Y;,t\Y n , 0) exp (-^k nx Y 2 )) S x (l £ fiY{ + l) S x (l £ / m y m + L) 

n=l / \ 1=1 / \ m=l x 



I! / dY n exp {~^k nx Y 2 )) S x (l £ f m Y m + L)) 

n=l / V m=l / 

and C z (t) are defined by 

2M \ / 2M \ f 2M 

]1 fdY>fdY n G(Y>,t\Y n ,0)exp(-^k nx Y 2 ) ) S y (L £ giY{)S y [L £ ^ 

n=M+l J \ l=M+l J \ m=M+l 



Cy(t) 



2Af \ / 2N 

n JdY n exp(-§-k nx Y 2 ) )S y [L £ g m Y n 

=AT+l J \ m=M+l 
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II JdY;fdY n G(Y>,t\Y n ,0)exp(-llk nx Yi) ) S z [l £ h{Y{ ) S z (L £ 9mK 

=2Af+l / V l=2M+l J \ m=2M+l 

II / dY n exp (-^rk nx Y%) )S Z [L £ h m Y m ) 

=2A/+1 / V m=2jV+l / 



cm 



On performing the integrations, one gets C x (t),C y (t) and C z (t) given in equations fl90|) to 

(E2D. 

APPENDIX D: THE EVALUATION OF S x (0), S y (0) AND 5,(0) 



The value of S x (0) can be found as follows. Defining the vectors f = (/i, fi-, ■■■■fj\r) an d 
a = (oi, a 3 , ....a 2 AA-i), £^(0) is given by 

at p2 
k. 



AT , 2 



n=l 



Using Equations (I4b[) and (To4|) . we get 

5,(0) = a.(Tf )- 1 /2 U 0o ( , u W( T J )-i/2 V ^( T ^-i/2 U ^ - 1 uW( T J°)-i/2. a t 

a. (Vfl .a' 



2a n _ 3 
^— ' n 2 7r 2 
n=odd 



8vr 2 ' 



Similarly, using Eq. (jJBJ) and (j5Sl) gives 



= Eir- 



(T^ 6 ) -1 / 2 "^ 6 ( u^ e '''(T^ 6 )~ 1 / 2 V* e (T^ e )~ 1//2 U ( ' !,e J xjM^'j 1 ^ ) -1 / 2 



(Vf)^ 1 



(Dl) 



:d2) 
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8vr 2 



and in a similar fashion 



9e\-l 



4tt 2 ' 



(D3) 

(D4) 
(D5) 

(D6) 
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